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We attempt to describe the stress distributions of granular packings using lattice-based layer- 
by-layer stochastic models that satisfy the constraints of force and torque balance and non-tensile 
forces at each site. The inherent asymmetry in the layer-by-layer approach appears to lead to 
an asymmetric force distribution, in disagreement with both experiments and general symmetry 
considerations. The vertical force component probability distribution is robust and in agreement 
with predictions of the scalar q model while the distribution of horizontal force components is 
qualitatively different and depends on the details of implementation. 



I. INTRODUCTION 



Gaining an understanding of the inhomogeneous stress distribution in a granular packing is important both be- 
cause of the insight it may yield into failure mechanisms and as a first step towards elucidating the dynamics of such 
systems p[-p^. Although qualitative aspects of the inhomogeneities have been known for some time ||,|,14-1^, quan- 
titative experiments have been performed only recently for cylindrical packings of glass beads |20[ , 2D arrangements of 
optical fibers |15|, and shear cells of glass spheres These experiments, together with numerical simulations [p2|-p4|, 
have yielded new insight into the nature of stress inhomogeneities. Various properties of the stress distributions are 
consistently observed both in two- [^,|2^,|2^ and three-dimensional ||^,|2^,^ systems. One such feature is that P{f), 
the probability of observing a force /, decays exponentially with / for forces much larger than the mean p|,p0|,p2[. 

A statistical model for a single component of stress in a packing based on a layer-by-layer approach appears to 
capture some features of the observed fluctuations In this model, the disorder in the system leading to the 

creation of force chains, whether from the inhomogeneities in the packing itself, variations in the size of the particles, 
or differences in material properties, is encapsulated in a set of random variables labeled qij which determine the 
fraction of the stress component is transferred from an element i of the packing to a neighboring element j. In ^ 
and 1^, these fractions are chosen randomly from probability distributions consistent with the constraint of force 
balance with the assumption that the g's at different sites are uncorrelated. Correlated q models have also been 
investigated |^,|2^. The layer-by-layer structure enables analytic progress on characterizing the force distributions. 
The predicted exponential decay in the probability distribution for large values of the stress component is in qualitative 
agreement with experimental and simulation results. 

The q model is scalar: it ignores the contributions that balancing the remaining stress components and torque 
may have in determining the weight redistribution and yields no information on their probability distributions. In- 
vestigation of correlations between fluctuations of shear and of compression effects, processes which are important 
in understanding failure, are not possible with this model. Moreover, though the q model successfully describes the 
stress fluctuations in simple geometries where the large-scale stresses are spatially constant, if applied to situations 
where stress varies over long scales, it predicts that these variations should obey a diffusion equation, in disagreement 
with experiment |^,|l^. These shortcomings have led several groups to examine stress fluctuations in models which 
incorporate vector forces [p7|-pl| . 

The generalization of the scalar q model to vector forces is an important step in the understanding of fluctuations 
about locally-averaged quantities calculated using continuum theories. Claudin et al. pd| ] have investigated the 
connection between the lattice-based q model and the "light-cone" continuum equations of Wittmer et al. by 
examining continuum equations with randomness. This connection is not trivial. In addition to subtleties encountered 
when one takes the continuum limit of stochastic models, Claudin et al. expose an important complication that arises 
when significant randomness is introduced into their continuum equations — the occurrence of tensile forces. They 
interpret this result, quite reasonably, as indicating that the granular material must rearrange. However, the materials 
should eventually reach a state at which the load should be supported and all the forces non-tensile. It is unclear 
how to describe this state using their approach. Generalization of the lattice-based models provides another means 
of approaching the issue. 
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Describing the system using random variables subject to the constraints of force balance at each site and no tensile 
forces provides a mechanism by which more realistic vector models may be based. Vector models in this spirit have 
been proposed by Eloy and Clement ||2^ and Socolar Common to the proposals is the propagation of forces 

downward in the packing in a layer-by-layer fashion starting at a load applied at its top. 

In this paper we present our attempts to construct a layer- by- layer vector force model. Wc find that serious 
fundamental problems arise from attempts to describe vector force fluctuations using simple generalizations of the 
q model. In particular, we find that it is very difficult to construct a stochastic model that leads to isotropic force 
fluctuations and satisfies the constraints of force and torque balance with non-tensile forces. By symmetry, isotropy 
of these fiuctuations is expected when a system is both prepared and compressed isotropically; moreover, Mueth et 
al. [20I demonstrate experimentally that subjecting a granular system to a uniform isotropic pressure leads to a stress 
fluctuation distribution which appears to be isotropic. The difficulty arises because of the inherent asymmetry in 
the formulation of a layer-by-layer vector model. At the individual site level, we see the possibility of the creation 
of large-magnitude horizontal output forces that are independent of the input forces and torque. The layer-by-layer 
structure does not allow for an intrinsic mechanism of removal, leading the most natural formulations of the model to 
have horizontal forces that appear to grow without bound as the system depth is increased. Incorporating cutoffs on 
the force magnitudes appears necessary to achieve any set of reasonable force redistributions. We find the probability 
distributions of vertical forces are insensitive to the cutoffs while the probability distribution of horizontal forces are 
dependent on their form. Thus, these cutoffs do not appear to present a solution to the underlying pathology of the 
model. 

The paper is organized as follows. Section || defines the model we investigate. We consider the redistribution of 
forces at a single site and discuss the difficulties that may be seen even at this level. Section III describes the process 
of redistributing forces through a lattice including our strategies for limiting the magnitudes of the forces generated by 
the redistribution algorithm. Section [V reports the results of the statistics of the forces obtained for forces obtained 
for various choices of the coefficient of friction and force cutoffs, illustrating the inherent asymmetry and cutoff 
dependence. Section ^ compares the results to experiment. Appendix ^ describes methods used to increase the 
efficiency of generating force redistributions. Appendix ^ discusses exactly solvable four-site lattice configurations to 
gain insight into strategies that could be used to construct isotropic models. 



II. MODEL 



A. Force Balance at a Site 

As in the scalar version of the q model we assume that the essential features of the disorder in a granular 

packing can be described using random variables. The choice of these variables is constrained by the requirements 
of satisfying force and torque balance as well as the requirement that the forces be non-tensile. We describe below 
the specific representation of random variables that we have employed. For the model to be considered useful, the 
probability distributions of vertical and horizontal force components should not be sensitive to the details of these 
choices. 

In our implementation of a layer-by-layer vector model we assume the topology of the packing is that of a modified 
regular two-dimensional triangular lattice of discs, as shown in figure |l|. Forces are introduced at the top of the lattice 
and are propagated downwards with "input" forces at a site arising from the two neighbors in the layer above and the 
resulting "output" forces being passed on to the two neighbors below. Sites on the same layer do not transfer force 
between one another. In an arbitrary N row by M column lattice, the jth site in the zth layer transmits its leftward 
output to the site j — 1 (j) on the i + I layer and its rightward to j (j + 1) for odd (even) i. 

We have chosen to parameterize the redistribution of forces at a site by randomly chosen contact angles and effective 
friction coefficients. Figure || shows a schematic of the forces acting on a site. The contact angles (pi,r determine the 
direction of the output normal forces Fi ^ to the left and right neighbors, respectively, and are chosen from the interval 
(0, -1). The effective friction coefficients rji ^ determine the direction and magnitude of the output tangential forces and 
are chosen from the interval [—1, 1]. The magnitude of the tangential force fd {d = l,r) varies from [0, j/irydi^dl] with fi 
being the coefficient of static friction. Positive values of rjd have tangential force fd contributing to the balancing of the 
vertical force component at a site in conjunction with its paired normal Fd-, while negative values have the tangential 
force in opposition to the paired normal. The tangential forces contribute to the torque TR at a site and are assumed 
to occur at the same distance R from its center of mass. The constraints imposed by force and torque balance and 
the non-tensile force requirement will place further restrictions on the allowed range of the random variables. 

In any stationary packing, the individual sites must satisfy force and torque balance: 
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F^" = —Fi{cosLpi — firji sinifi) + Fr{cos(pr — ^r],. sinipr), (la) 
fy" = Fiisimpi - ^rii sin ipi) + Fi{sm(pr - fn]r sin ipr), (lb) 
r™ = fiirjiFi - r,rFr), (Ic) 

where the total input force components and torque are given by the sum of the normal and frictional forces from 
neighboring sites in the layer above: 

= (Fr), + (Fj.")^, (2a) 
i^r = (Fi"), + (Fr"),' (2b) 

r" = rj" + r™. (2c) 

Only the total input force and torque enter into eqs. (|l|) as these are fixed values arising from the propagation of 
forces in the previous layer. Because of torque balance, eq. (px|), the 77's are not independent — r]i may be written in 
terms of rjr (or vice versa). The number of independently chosen random variables is reduced to three. 
Solving for the output normal forces yields: 

Fi = [F™ (sin ipr + urir (cos (pi + cos ipr)) - Fy"- (cos ipr + (sin tpi - sin Lpr)) 

+r" (cos {lpi + ipr) - fn]r sin {ipi + (fir))] , (3a) 
Fr^i [-F:^ sin - F;" cos + r"'] , (3b) 

where 

'i' = -^77r [1 + cos {lpi + ipr)] - sin {ipi + ip^) ■ (4) 

The asymmetry between the forms of Fi and Fr is due to the choice of designating rjr as a random statistical variable. 
The force redistribution is considered valid if for the values of ipi^r and rjr chosen, 

Fi^r > (Non— tensile constraint), (5a) 
^ ^ (Newtonian friction). (5b) 

Figure ^ shows the valid region of contact angle-effective friction coefficient configuration space for a vertical input 
force with positive and negative values of rjr- The formalism for choosing r]i as a random variable instead of rjr is 
similar. 



B. Difficulties Arising at a Site 

Even at the single-site level, we may encounter difficulties in the redistribution and propagation of forces. One 
possibility is that the input forces and torque may be such that no valid redistribution can occur. Although the force 
and torque balance equations always yield solutions for the inputs in any randomly chosen configuration of angles 
and friction, the additional non-tensile and frictional constraints may severely limit the number that may be realized. 
Typically, we see that the configuration space is significantly reduced if the net input force is largely horizontal or if 
the input torque is large relative to the input normal and frictional forces. 

A more serious difficulty arises because large magnitude output forces can be generated irrespective of the magnitude 
of the input force. To illustrate this problem, we write the force balance at a site, given by eqs. (Q), schematically as 

F-=M(^,,^r,?yr)F°:;^^,l, (6) 

which has the solution 

F°:;^^ai = M-^F-. (7) 

The factor from cq. (jj) is the determinant of the matrix M: 

* = det M. (8) 

As Fi^r oc 4*^^, large magnitude normal forces will be generated whenever ^' approaches zero. ^ depends only on the 
choice of random variables and is independent of the input force. 
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Exactly at the ^ = boundary, the output forces sum vectorially to zero — no force balance is possible. This 
boundary occurs at values of ipi^r and ?7r < satisfying 

The boundary can be clearly seen in figure ^(a). Although the boundary itself does not yield solutions, regions of 
configuration space exist near it that do yield valid solutions for which j^*! is very small. In these regions the output 
forces almost cancel, requiring large magnitude normal forces in order to satisfy any equilibrium condition for non-zero 
input values. Because the input force is designated as originating from the neighboring sites from the layer above, 
an asymmetry exists between vertical and horizontal components. The vertical components of the output forces are 
bounded in magnitude by the fixed vertical component of the input while the horizontal components of the leftward 
and rightward outputs acting in opposition to each other are essentially unbounded. This effect can easily be seen 
in the non-frictional case for a purely vertical input force: as shown in figure ^ as the angle of both output forces 
approach the horizontal, the magnitude of those forces must increase so that their vertical components will balance 
the input. The frictional case is similar, though the addition of tangential forces complicates the picture slightly. 

The vector model proposed by Eloy and Clement [ p7| also suffers from this pathology although the mechanism may 
not be as obvious due to their choice of parameterization. Socolar's model Eq] requires the forces to lie in a 45° cone 
about the vertical and therefore only considers non- negative values of friction (as interpreted by our model). As a 
result, it does not generate these large forces but at the cost of imposing severe limits on the magnitude of horizontal 
components. We discuss this issue in more detail below. 



III. METHODS 



A. Force Redistribution through a Lattice 

Our algorithm for generating large lattices redistributes the forces at individual sites starting from a normally 
distributed vertical load applied on the top layer of sites and proceeding downward into the lattice layer by layer. All 
sites on a layer are processed before proceeding to the next layer down. Contact angles ipi^r and an effective friction 
coefficient, either rji or ry,., at a site are randomly chosen within their respective intervals using a uniform deviate and 
a test is made to determine whether the resulting force redistribution satisfies the necessary non-tensile and frictional 
constraints. If the constraints are not met, new sets of random variables for the site are chosen until all requirements 
are satisfied. 

Failure in lattice generation occurs when the input forces and torque at a site cannot be be redistributed within 
a reasonable sampling of the contact angles-effective friction coefficient configuration space. Reasonable has been 
defined as a sampling of 25,000 uniformly distributed points in the space. A larger sampling size did not increase the 
lattice yield significantly. Furthermore, some input force configurations do not have any valid redistributions. In our 
simulation, if no valid redistribution for a site is found, lattice generation is terminated and restarted with another 
random number seed. While encountered mainly when friction has been applied, we find that non-frictional cases 
may also suffer from this behavior. The rate of failure increases with lattice size and coefficient of static friction /x. 
However, techniques described in appendix can be used to increase the yield. 



B. Implementation of Limits on Force Magnitude 

As discussed above, large horizontal forces can be generated at a single site. However, we consider the possibility 
that the lattice structure may be self- limiting so that the effect of large magnitude horizontal forces will be restricted 
to small neighborhoods surrounding the generating site. One possible mechanism would be the cancellation occurring 
at a site with left and right input neighbors both contributing this type of force. However, as shown below, we find 
that the generation of a cancellation pair is unlikely enough that these large forces build up and propagate in our 
lattices. 

Having failed to identify an intrinsic means to limit the production and transmission of large magnitude horizontal 
forces, we impose various cutoff schemes to attempt to generate sets of realistic force distributions. The choice of 
one form of implementation over another is somewhat arbitrary. Consequently, several cutoff schemes have been 
implemented to observe the influence they exert on the resulting force distributions. 

The first cutoff scheme we implement is the simplest: the magnitude of each normal force at a site may not exceed a 
specified value. The second scheme is to limit the allowed angles to prevent exploration of the = boundary. This 
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method is similar to the Eloy and Clement model ||2j] which has fixed angles. By restricting the range of available 
angles, this type of limit serves to reduce the volume of configuration space available for redistribution and can exclude 
regions that form large magnitude horizontal forces. Our third scheme is a soft cutoff scheme based on the assumption 
of contact energies between sites following a Boltzmann-like distribution. This choice is clearly arbitrary as the system 
is not thermal. For simplicity, the lattice is modeled as a layer of spheres whose centers lie on the same plane (contact 
topology is therefore equivalent to discs) and elastic theory is used to calculate the energy U within the contacts 
between these spheres: 

U = A{a,E,R)F^, (10) 

where F is the normal force between the sites in contact and ^ is a function of the material properties (Poisson's 
ratio a and Young's modulus E) and radius R of the spheres. We assume that the probability that a contact has an 
energy U follows an exponential distribution 

P(C/) = -^e-^/^°, (11) 
where Uq is is an arbitrarily assigned average contact energy. 

C. Generation of Datasets 

For each coefficient of friction /i and force cutoff configuration a set of 1000 horizontally periodic lattices of 100 rows 
by 100 columns with an applied load of 1000 N (the unit is arbitrarily applied for the benefit of the applied limits) is 
generated and averaged over to yield the probability distributions of vertical and horizontal force components. The 
load is distributed on the top layer via a normal distribution centered about the average force of 10 N with a standard 
deviation of 5 N. Normalization of both the vertical and horizontal force components is performed relative to the 
average vertical force (10 N) imposed on a site by the load on the system. The lattice size is chosen to minimize 
computational time as the failure rate (number of lattices which failed to run to completion divided by number of 
lattices started) in lattice generation grew significantly with increased vertical size despite the measures taken to 
increase overall yield described in appendix However, as convergence of the distributions proves to be fairly rapid, 
larger lattices are unnecessary. 

Data sets with static coefficients of friction fi = 0, fi = 0.1, and fi = 0.2 under various cutoff schemes are used. 
Configurations for /x = and /i = 0.2 are also generated without an applied force cutoff in order to demonstrate 
the necessity of cutoff implementation. In setting the sharp force cutoff, upper limits of 50 N and 100 N are set on 
the normal forces. Angle cutoffs are implemented for a range about 60° to simulate a triangular packing. For the 
soft cutoff, eq. (pi]), we use the values based on the physical properties of soda hme-sihca float glass (cr = 0.23, 
E ^ 7.2x 10^° Pa) assuming a radius R of 1.75 mm, yielding A{a,E,R) = 1.5 x lO^'^J • N^i for eq. (0). Total 
input energies of 1 and 5 J are considered, leading to Uq = 5 x 10~^ J and Uq ~ 2.5 x 10^* J, respectively. 

IV. RESULTS 

We find the probability distribution P{v) of normalized vertical forces v shows remarkable robustness and appears 
to be independent of the coefficient of friction ^ and of choices of large force cutoff. In contrast, the probability 
distribution P{h) of normalized horizontal forces h exhibits changes in functional form with variation in both ji and 
with cutoff choice. Horizontal forces, in general, are of larger magnitude than vertical forces. 

Convergence of the probability distributions of normalized vertical force component v is fairly rapid, on the order of 
10 rows. P(v) versus depth in lattice is shown in figure |(a) for a non-frictional, sharp cutoff limit configuration. Similar 
results for P{v) are seen for all generated configurations with little variation in functional form; the distributions at 
row 100 are shown in figure |^(b) with the symbol key for the various configurations shown in table ^ An exponential 
tail for P{v) is seen for larger values of v and a "dip" in P{v) is seen for small values of v. The observed P{v) is 
very similar to that obtained from the scalar q model with uniform q distribution for an iV = 2 (two-dimensional) 
system [|[ 

P{v) = Av'^e-^\ (12) 

The distributions of q values describing the redistribution of the vertical forces, shown in figure ^, are nearly uniform 
with a slight increase in probability near values of and 1 as compared to 0.5 for fi = and a decrease when /x > 0. 
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In addition, configurations with /i > in the vector model do allow for q values outside the range [0, 1] due to mostly 
horizontal normal forces with a corresponding 77 < 0. The net vertical component of the input force at a site is still 
kept positive. For the fJ, configurations examined, roughly 5-20% of the q values are found to lie outside the 
range [0, 1] with a smaller fraction found in more restrictive force cutoff configurations and a larger fraction for less 
restrictive cutoffs. Non-frictional (/i = 0) configurations all have q values lying in the range [0, 1]. 

In contrast to the vertical force distribution P{v), the probability distribution P{h) of normalized horizontal force 
component h shows great variation in form with changes in the imposed force limit. We first examine the P{h) 
distributions for lattice sets with no imposed force cutoff and = and fj, = 0.2, shown in figure 0. We see the 
distributions spreading toward larger values of h as we progress deeper into the lattice while P{v) remains bounded 
and robust. When /i = 0, P{h) converges, but at values of the horizontal force much greater than the vertical force, as 
seen in |^(a) . We find that this convergence is mainly to a decrease in the number of valid solutions of force and torque 
balance, eqs. (|l|), that satisfy the non-tensile and frictional constraints given by eqs. as the ratio of input force 
components F^"/F™ grows larger. Although this process appears to provide an intrinsic limit on forces, it results in 
physically unreasonable values of force and fails altogether to limit the magnitude of horizontal forces when fi > 0. 
The addition of friction greatly enhances the anisotropy — the generation of large horizontal forces in configurations 
without friction only occurs at nearly horizontal values of fi.r, while the addition of friction allows for the generation 
to occur for a larger region of angles. More precisely, the 5* = boundary for /i = occurs only at 

ipi + ifj. ^ 0, (13) 

while it varies for /i = 0.2 according to 

< ipi + Lpr < 22.6°. (14) 

The rapid growth in the magnitude of h values prevents the generation of large frictional lattices without imposed 
force cutoffs to take place in a reasonable amount of time. When a sharp cutoff in the normal force is imposed, the 
probability distribution P(h), as seen in figure ^(a) cuts off abruptly, which is unlikely to be realized in a physical 
packing. The use of angle cutoffs was unsatisfactory as lattice generation for this scheme was consistently terminated 
due to angles being driven to the cutoff values due to our choice of parameterization and random value selection. 
Configurations obtained using the soft cutoff from eq. ([ll|), shown in figure ||(b), exhibit differing functional forms 
for the probability distribution as the cutoff is changed with broadening occurring at the tail with increased input 
energy. For small values of h, the distribution appears to be unaffected by changes in the choice of cutoff. Increasing 
H serves to broaden the distribution P{h). As force components are normalized by the same factor, it can be readily 
seen that the scale of horizontal forces is larger than the vertical as shown in the inset of ||(a). A grayscale plot of the 
forces on a representative sample lattice with fi = 0.2 and a sharp cutoff at 100 N is shown in figure ^. It is obvious 
that the vertical and horizontal force fiuctuations differ qualitatively. 

Thus, although the imposition of force limits results in horizontal force distributions which do not diverge as the 
depth is increased, these limits serve only to mask the divergent behavior of the model: the distribution of forces 
expands to fill the space allowed and the distribution of horizontal forces exhibits a strong dependence on the choice 
of cutoff scheme and value. 



V. DISCUSSION 

The layer-by-layer vector model that we have investigated to model forces in granular packings yields vertical and 
horizontal force probability distributions which are of different functional forms and scales, in contrast with the recent 
measurements by Mueth et al. |2^]. Moreover, while the distribution of vertical forces is robust, the distribution of 
horizontal forces depends on the details of our implementation. 

The asymmetry seen in the force distributions is a reflection of the vertical-horizontal asymmetry inherent in any 
uni-directional layer- by-layer model. At the level of individual elements, we see that large-magnitude horizontal forces 
are generated. As these forces accumulate, the resulting probability distribution of horizontal force components is 
naturally skewed toward larger values. In contrast, large vertical forces are not generated. The imposition of limits 
on the vertical forces by the applied load and the non-tensile force constraint is sufficient to limit the magnitude of 
that force component. The horizontal component has no similar constraints except those arbitrarily imposed by our 
cutoff schemes. 

This behavior is apparent in other implementations of layer-by-layer models despite differing choices in parame- 
terization and configuration. Eloy and Clement [ p7| model the packing by assuming a mono-disperse 2D array of 
hard cylinders arranged in a triangular lattice. The angle of contact between cylinders in neighboring layers is fixed 
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at slightly less than 60° measured with respect to the horizontal. The redistribution of forces at a cylinder is pa- 
rameterized by the coefficient of friction fj, and the difference in value of the horizontal force components transfered 
to the neighboring sites in the row below, p. They clearly note that for certain values of /i, valid values of p are 
unreasonably large in magnitude and an arbitrary cutoff, restricting the parameter space to be far away from this 
divergent region, is imposed. Socolar's a model |28| represents the packing with a lattice of square cells with net 
normal forces, couples (torques), and tangential forces represented on each edge. The redistribution of forces at a 
cell is parameterized by a triplet of values {ao, ai, 02) which are used to couple the torque and force balance. The 
distribution of net forces at every site is restricted to be in a downward cone opening at 45° by considering only 
frictional forces acting in conjunction (the same direction vertically) with normal forces. This restriction guarantees 
successful lattice generation but renders the force fluctuations anisotropic by construction — the magnitude of the 
horizontal force components is artificially limited to less than that of the vertical. We believe that widening the cone 
by including frictional forces in opposition to normal forces would lead to the divergent behavior of horizontal force 
components. 

Transmission of information regarding the production of large-magnitude forces upwards into the already processed 
portion of the lattice is a mechanism missing in layer-by-layer models which could serve as means of suppressing large 
magnitude forces. To illustrate this point, figure nu shows the breakup of a large horizontal force into upward- and 
downward-traveling components . In addition to reducing the magnitude of a single force which would have otherwise 
persisted, this process yields an upward force component that may instigate a rearrangement which then causes a 
new redistribution of forces in the layers above and allow the site generating the force to "relax" to a less stressful 
configuration or generate a loop which isolates the force chain. One must consider some sort of feedback mechanism 
to enable corrective rearrangements and redistribution of forces to be made. 

Implementation of a feedback mechanism within the context of a layer-by-layer method is not a trivial matter. While 
a pseudo-feedback mechanism (see appendix ^ is used in the generation of the sets of lattices analyzed here and 
serves the purpose of increasing the likelihood of finding an allowable configuration, it still maintains the downward 
propagation of actual force information. 

In summary, the intrinsic asymmetry of layer-by-layer vector models results in asymmetries in force distributions 
that are not seen in experiment. This failure of the model indicates that constructing a statistical characteriza- 
tion of vector force transmission in granular media requires a better understanding of the vectorial nature of force 
redistribution, taking into account explicitly the symmetry properties of the medium. 
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APPENDIX A: PSEUDO-FEEDBACK AND INCREASED LATTICE YIELD 

The generation of a sufficiently large set of lattices is essential in performing any reasonable statistical analysis of 
the force distribution in the model. In the q model, valid redistribution of force at a site is guaranteed because the 
constraints on the choices of g's are independent of the force acting at a given site. In the vector model proposed here, 
as well as in and |2^] , the set of random statistical variables allowed at a site is dependent on the force acting upon 
it, as can be seen by eqs. (||) and (|^) for calculating the output normal forces and checking the non-tensile constraint, 
respectively. We choose to select the set of random statistical variables for a site, test to see if the constraints are 
met, and discard and reselect a new set of variables if the constraints are not satisfied. Although seemingly inefficient, 
this method is far simpler than encapsulating the non-tensile and friction constraints within the choice of random 
variables. The lack of a guaranteed redistribution configuration is the root cause of the failure in lattice generation 
in the vector models we examined. 

We may reduce the rate of incidence of failure if we can ensure that the inputs at the sites in the layer below 
that currently being redistributed will be able to support valid force redistributions themselves. We do this by 
incorporating a simple check in the selection process for the statistical variables (fi^r and rji^r for neighboring sites in 
the same layer. If a valid redistribution is not possible, then all sites on the row are subjected to a new redistribution. 
We may extend this process to an arbitrary number of resulting layers although this may significantly increase the 
memory requirements for computation in addition to complicating the layer-by-layer algorithm. The essence of this 
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requirement is to capture the flavor of feedback and rearrangement, albeit in an imprecise manner. However, we have 
no guarantees that a valid configuration will be possible. 

This pseudo-feedback (PFB) addition was set to a depth of two layers in the process of generating the lattices used 
in this paper. Typically, most difficulties in redistribution may be resolved with a one-layer PFB as this provides an 
immediate check for input forces leading to a valid redistribution. However, these valid redistributions may lead to 
an invalid input force on the next layer down, a less likely but still significant cause of failure. The implementation 
of a second layer for the PFB resolves this issue and allows for the vast majority of lattices for our configurations of 
friction and force cutoff to run to completion. Using PFB beyond two layers increases the yield but not significantly. 



APPENDIX B: CONSIDERATIONS RELEVANT FOR SYMMETRIC LATTICES: FOUR-SITE LATTICES 

We have examined exactly solvable lattice configurations to study the interaction between force-balanced sites 
to gain insight into how large-magnitude forces may persist in the lattice and what methods beyond externally 
imposed force limits exist to prevent or remove them. Although they do not offer any statistical information, these 
configurations offer further insight into the nature of the generation and propagation of large magnitude horizontal 
forces by allowing greater control over the parameters governing the forces in a packing than is found in a stochastic 
layer-by-layer approach. 

We find that four-site lattices represent the smallest meaningful unit of study. The two configurations which 
maintain the coordination number four are the two-layer triangular lattice and the quadrilateral shown in figures ^ 
and respectively. We examine both horizontally periodic and fixed-wall cases. 



1. Two-Layer Triangular Lattice 

Placing the sites in the familiar triangular lattice shown in figure ^ vertical forces are applied to each site. 
Because of the symmetry of the system, two angles and a single internal friction coefficient are sufficient to specify 
the redistribution of forces for the internal contacts. Assuming that the known angles are ifis and (^23 and the known 
friction coefRcient is given by 77 = 7713, the normal forces for the periodic case are 

-^"13 = cos (p23 (Sin((/3i3 + (P23)+ m (1 + C0s((/?i3 V323)))~\ (Bla) 

F23 = Fi3 (cosv5i3 - fj,T] (sin (^13 - sin (^23 ) ) / cos (/923 , (Bib) 
F24 = Fi3 (cos ipi3 - ^irj (sin (^13 - sin (^24)) / cos (^24, (Blc) 
Fii = Fi3 (cos Lpi3 - ^irj (sin (^13 - sin ipu)) / cos (^14, (Bid) 

where (/?i4 and <^24 are defined implicitly by 

FT_ _ cosy24 (sin(y)i3 ^23) + M'? (1 + cos(y)i3 + ^23))) .^2 ] 

^2 cos ipx3 Sin((p23 + <^24) + m (cos ip23 (1 + COs((^23 + V24)) - ) 

FT_ ^ COS ip23 (sin(y)i3 + (^14) + m?? (1 + cos(tpi3 H- 

F3™ cos (/7i4 (sin(.^i3 + (^23) + M?? (1 -I- C0s((pi3 -I- ip23))) ' 



Assuming uniform vertical inputs, the periodic case admits solutions near the single-site 4* = boundary with the 
balancing of external forces remaining perfectly satisfied while the magnitudes of the internal normal and frictional 
forces grow unbounded. A similar result exists in the non-uniform case once external force and torque balance is 
taken into account. 

The fixed horizontal boundaries case is derived by disconnecting sites 1 and 4 and fixing the direction of the normal 
forces to be normal to the wall. As the magnitude of the sum of the horizontal components of each normal and 
frictional force for each contact must be equal, any restriction placed on the at the wall (ie. loading) will place 
limits on the internal redistribution. By removing the asymmetry of loading in this configuration, we can prevent the 
formation of the large magnitude forces. However, we must be able to specify all the forces along the boundary of the 
system. If we are only given the vertical loading, then any restrictions placed on the horizontal force component will 
be arbitrary. 

This result may indicate that successful implementation of probabilistic vector models of force fluctuations requires 
better understanding of the role of the boundaries. This question is also crucial to understanding whether the equations 



underlying these force distributions are elliptic or hyperbolic |11 
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2. Quadrilateral Lattice 



The quadrilateral configuration shown in figure 12 offers the advantage of having a "rotational" symmetry lacking 
in both the single-site and triangular lattice configurations. We have implemented an explicit loop to observe its effect 
on the internal force redistributions and hopefully allow us insight into the role of this symmetry in the system. 

We apply vertical external forces to each site. Periodic boundary conditions are created by pairing the horizontal 
inputs of sites on the same vertical "level" so that they are equal in magnitude and have friction coefficients acting on 
each site of the pair equally and oppositely. Although the equations of force balance for the system are easily derived, 
valid configurations are more readily determined numerically. Unbounded horizontal force solutions are supported — 
the redistribution pair consists of the periodic horizontal and the internal force connecting same-level sites. The 
internal loop prevents these forces from being passed between vertical levels due to the direction of the tangential 
frictional forces within the loop being prevented from aligning adversely in neighboring loop sites. 

The addition of fixed horizontal boundaries leads to results similar to the triangular lattice configuration. We see 
that the introduction of internal symmetry to the system is not enough — we are still subject to the effects of an 
externally imposed asymmetry. However, we do gain the ability to isolate the effects, in this case over a layer, as 
opposed to its propagation throughout the system. 



[1] C. Liu et al, Science 269, 513 (1995). 

[2] S. N. Coppersmith et al, Phys. Rev. E 53, 4673 (1996). 

[3] P. Dantu, Ann. Ponts Chauss. IV, 193 (1967). 

[4] A. Drescher and G. de Josselin de Jong, Jnl. Mech. Phys. Solids 20, 337 (1972). 

[5] T. Travers et al, Europhys. Lett. 4, 329 (1987). 

[6] S. F. Edwards and C. C. Mounfield, Physica A 226, 1 (1996). 

[7] C. C. Mounfield and S. F. Edwards, Physica A 226, 12 (1996). 

[8] S. F. Edwards and C. C. Mounfield, Physica A 226, 25 (1996). 

[9] J. P. Wittmer, P. Claudin, M. E. Gates, and J.-P. Bouchaud, Nature 382, 336 (1996). 
[10] J. P. Wittmer, M. E. Gates, and P. Glaudin, J. Phys. I France 7, 39 (1997). 
[11] P. Glaudin, J.-P. Bouchaud, M. E. Gates, and J. P. Wittmer, Phys. Rev. E 57, 4441 (1998). 
[12] J. Hemmingsson, Physica A 230, 329 (1996). 

[13] J. Hemmingsson, H. J. Herrmann, and S. Roux, J. Phys. I France 7, 291 (1997). 
[14] M. Ammi, D. Bideau, and J. P. Troadec, J. Phys. D 20, 424 (1987). 

[15] G. W. Baxter, in Powders and Grains 97, edited by R. P. Behringer and J. T. Jenkins (Balkena, Rotterdam, 1997), pp. 
345-348. 

[16] P. Dantu, in Proceedings of the 4th International Conference on Soil Mechanics and Foundation Engineering, London, 1957 

(Butterworths, London, 1958), pp. 144-148. 
[17] T. Wakabayashi, in Proceedings of the 9th Japan National Congress for Applied Mechanics, Tokyo, 1959 (Science Gouncil 

of Japan, Tokyo, 1960), pp. 133-140. 
[18] T. Travers et al, J. Phys. (France) 49, 939 (1988). 

[19] D. Howell and R. P. Behringer j_in Powders and Crains 97, edited by R. P. Behringer and J. T. Jenkins (Balkena, Rotterdam, 

1997), pp. 337-340. See Ref. ||. 
[20] D. M. Mueth, H. M. Jaeger, and S. R. Nagel, Phys. Rev. E 57, 3164 (1998). 
[21] B. Miller, G. O'Hern, and R. P. Behringer, Phys. Rev. Lett. 77, 3110 (1996). 
[22] F. Radjai, M. Jean, J. J. Moreau, and S. Roux, Phys. Rev. Lett. 77, 274 (1996). 
[23] G. Thornton and S. J. Antony, submitted to Phil Trans Roy Soc A (1998). 
[24] R. Glelland (unpublished). 

[25] P. Glaudin and J.-P. Bouchaud, Phys. Rev. Lett. 78, 231 (1997). 

[26] M. Nicodemi, Phys. Rev. Lett. 80, 1340 (1998). 

[27] G. Eloy and G. Glement, J. Phys. I France 7, 1541 (1997). 

[28] J. E. S. Socolar, Phys. Rev. E 57, 3204 (1998). 

[29] E. B. Pitman, Phys. Rev. E 57, 3170 (1998). 

[30] V. M. Kenkre, J. E. Scott, E. A. Pease, and A. J. Hurd, Phys. Rev. E 57, 5841 (1998). 
[31] J. E. Scott, V. M. Kenkre, and A. J. Hurd, Phys. Rev. E 57, 5850 (1998). 

[32] L. D. Landau and E. M. Lifshitz, Theory of Elasticity, Vol. 7 of Course of Theoretical Physics, 3rd ed. (Butterworth- 
Heinemann, Oxford, 1986). 



9 



Row 




FIG. 1. A triangular lattice is used to approximate a 2-D granular packing in order to assign depth in site and neighboring 
sites. Lines connecting sites indicate that a transfer of force occurs between them. Sites on the same layer are not connected. 
Site {i,j) is labeled along with its neighboring sites. 




FIG. 2. Schematic of forces at a site for the vector model. Input forces are from the neighboring sites above and output 
forces are from the sites below. The frictional forces have magnitude equal to fi^r = f^Vi.rFi^r, where Fi^r is the normal force, 
and are shown in the +7] direction. The angles ipi^,. indicate the contact angle between the site and the output neighbors. The 
parameters ipi^r and rii,r are chosen randomly, consistent with the constraints of force and torque balance, non-tensile forces, 
and |r7i,r| < 1- 
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(a) (b) 

FIG. 3. Slices of configuration space of contact angles (ySj.r and effective friction coefficient rjr. The maximum coefficient of 
static friction fi has been set to 1 to exaggerate the features. A sample of 10,000 points were randomly chosen with uniform 
probability in the interval [0, ^] x [0, Points that satisfied the non-tensile and frictional constraints, eqs. are marked 
with a +. Boundaries of the valid region are labeled. The $ = boundary where large horizontal forces are generated is found 
for rjr < but not for ry^ > 0. 



9=45° 



9= 26.6* 




FIG. 4. Even for a site with a purely vertical input force in the absence of friction, the output forces must increase significantly 
in magnitude as both output angles approach the horizontal. Because the horizontal force components in the outputs are in 
opposition to each other, they may grow without bound while the vertical components are limited in magnitude by the input. 
We see this growth from a small magnitude in (a) to a large magnitude in (c) for the same vertical input force. Values of 'I' 
from eq. (0) for each case are 4'45<' = —1, *26.6° = —0.8, and *i4° = —0.47. 
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(a) (b) 

FIG. 5. (a) Probability distribution P{v) of normalized vertical force v at increasing depths of a non-frictional (/i = 0) lattice. 
Force limit has been arbitrarily set to 100 N (unnormalized) for normal forces. Convergence occurs fairly rapidly (within ~ 10 
layers). Similar results are seen for other values of fj, and force cutoff configurations, (b) Probability distributions P{v) for 
normalized vertical force v at depth 100 for various values of /x and force cutoff configurations. The symbols are defined in 
table |. Functional form of the distributions is invariant with respect to configuration. The force distributions are very similar 
to that of the scalar q model with N — 2, shown as the solid line, which is appropriate for this geometry M. 
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FIG. 6. Histogram of qi values describing the redistribution of the vertical component of force for various values of and 
force cutoff configurations at a depth of 100 rows, denoting the fraction of weight supported by the leftward neighbor of a site. 
Symbols used are the same as for figure ^b) (defined in table |) except only the larger, less restrictive cutoff values are used 
(Sharp Cutoff at 100 N and Soft Cutoff with 5 J). Bin size, Ag, is 0.04 and the result is placed at the right edge of the bin. 
The distribution is nearly uniform except for a slight increase in probability near values of and 1 as compared to 0.5. for 
non-frictional configurations and a slight decrease when friction is present. Because of the admission of friction in the vector 
model, q values outside of the range [0, 1] are possible as frictional forces may cause a contact to have an upward net force. 
Roughly 5-20% of g's in /i 7^ configurations are found to lie outside the range [0, 1] with a smaller number found for more 
restrictive force cutoff configurations and a larger number for less restrictive cutoffs. Non-frictional (/i = 0) configurations all 
have q values lying in the range [0, 1]. 
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(a) (b) 

FIG. 7. Results for lattices without imposed force cutoff. We see that the probability distribution P{h) for normalized 
horizontal forces h broadens and shifts toward larger values of h with depth while P{v) for normalized vertical forces v remains 
unchanged. The distribution P{h) in the non-frictional {fi — 0) case appears to converge after about 100 layers, reaching 
upper limit of (h) / (v) ~ 5.5, as shown in the inset. The distribution does not widen any further — for large Fl^" /Fy" ratios 
at individual sites, no valid solutions of force and torque balance, eqs. (|l]), under the constraint of non-tensile forces, eq. (pa]), 
exist. Although this appears to be an intrinsic limit, it results in physically unreasonable values of force and is far less effective 
when /i > 0. Results shown in (b) demonstrates that the addition of friction greatly enhances the anisotropy. Lattice generation 
with friction present was limited to 5 rows due to an increased rate of failure. Deeper lattices are possible; however, the number 
of points sampled in configuration space increases significantly. 
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(a) (b) 

FIG. 8. (a) Probability distribution P{h) of normalized horizontal force h at increasing depths for a non-frictional (p = 0) 
lattice. The normalization factor used is the same as for the vertical component of force. Force limit has been arbitrarily set 
to 100 N (unnormalized) for normal forces and the distribution terminates abruptly near h = IQ (100 N). The distribution has 
not converged and appears to be widening, indicating the average horizontal force is increasing with depth. The functional 
form is unlike the vertical and is on a larger scale as can be seen in the inset, (b) Probability distribution P[h) at depth 100 for 
various values of /i and force cutoff configurations. Symbols used are the same as for figure |^(b) (defined in table |). Functional 
form varies with the imposed limits. While P{h) curves for the sharp cutoff terminates at /i = 5 and h = 10, the soft cutoff 
Boltzmann-liko limits do not. 
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(a) q Model (b) Vector Model: Vortical Force Com- (c) Vector Model: Horizontal Force 

ponent Component 

FIG. 9. Grayscale plot of representative 100 by 100 lattices for the q model with uniform q distribution and the vector 
model with ^ = 0.2 and a sharp cutoff of 100 N for normal forces. A load of 1000 N has been applied at the top of the 
lattices. The darkness or brightness of a site corresponds to the magnitude of the normalized force component. Normalized 
force components with magnitudes greater than 5 have been clipped at 5. Both the vertical and horizontal force components 
are normalized by the same value. The qualitative difference between the two components for the vector model are readily 
apparent; the horizontal fluctuations are of a much larger scale and form "V" -shaped chains reminiscent of light cones while 
the vertical fluctuations resemble those found in the q model. 



Without Feedback 



With Feedback 




FIG. 10. Schematic of feedback arising from a large horizontal force. The left-hand side shows the force propagation as 
implemented in layer-by-layer models. The right-hand side demonstrates how the force could be broken into upward- and 
downward-traveling components. 
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FIG. 11. Schematic of the four-site triangular lattice. External vertical forces have been applied to each site. The configu- 
ration is made periodic by connecting sites 1 and 4 as shown. Indexing of the contact angles and effective friction coefficients 
refer to the sites in contact (ie. 9513 is the contact angle between sites 1 and 3). As angles are mesisured with respect to the 
horizontal, <^ij = ifiji. 




FIG. 12. Schematic of the four-site quadrilateral lattice. External vertical forces have been applied to each site. The 
configuration is made periodic by connecting sites on the same vertical layer as shown (site 1 to site 2, site 3 to site 4). Contact 
angles are measured by their deviation from a square (ie. (^12 is the deviation from the horizontal and (^13 is the deviation 
from the vertical and ifij = —<fiji)- 
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TABLE I. Symbol key for figures ^(b), and ^(b) which are plots consisting of various configurations of coefficient of static 
friction ^ values and force cutoff schemes. 
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